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GEOSTAR-H 

A GEOPOTENTIAL AND STATION POSITION RECOVERY SYSTEM 

C. E. Velez 
G, P. Brodsky 

Program Systems Branch 
Mission And Trajectory Analysis Division 


ABSTRACT 

The GEOSTAR-H multiple arc geopotential co- 
efficient and station position recovery system 
is described. This system was developed using 
the GEOSTAR-I system as documented mNASA 
X-533-69-544, as a base, and is the second in a 
series of geodetic programs under current de- 
velopment and use by the Mission and Trajec- 
tory Analysis Division. 
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GEOSTAR-H 


A GEOPOTENTIAL AND STATION POSITION RECOVERY SYSTEM 


I. INTRODUCTION 

The GEOSTAR-II system is a multiple arc, multiple satellite geopotential co- 
efficient and station position recovery system. The basic features of this system 
are: 

® The capability to form the normal equations for up to 250 geodetic and 
tracking model parameters per arc 

• The capability to utilize a-priori {information m the formation of a 
multiple arc least squares solution for up to 500 parameters 

• The use of advanced numerical integration techniques 

• The availability of a "stand-alone" orbit generator program which can 
be used to measure the performance of the mtegration process used for 
the equations of motion and variational equations . 

The basic operational aspects, mathematical techniques and system design 
characteristics of the GEOSTAR-II system are based on the GEOSTAR-I sys- 
tem as documentedmNASAX-553-6 9-544. In particular, all the current capabili- 
ties of the GEOSTAR-I system are available for use in the new system; however, 
because of additional core stor age requirements of the n system, it is not recom- 
mended for use m lieu of the I system unless one or more of the new features 
are required, such as the extension to 250 parameters per arc. 

The GEOSTAR-II system is written in 360 Fortran IV and is currently operational 
on the 360/95. The operating instructions are given in reference 4. 

H. SYSTEM DESCRIPTION 

The primary module of the GEOSTAR-n system is a modified version of GEO- 
STAR-I ODP program which has been extended to handle up to 250 model 
parameters per arc. Other modifications allow the use a-priori estimates and 
statistics m the multiple arc mode and the use of a more rapid mtegration start- 
ing process based on an iterative scheme using closed form Newtonian type for- 
mulas. Other system programs include MERGE, SOLVE, EIGENVALUE and 
ORBGEN. 
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The GEOSTAR-II system can be used to form the normal equations required to 
estimate up to 500 orbital and geodetic parameters using any number of tracking 
data arcs. This can be accomplished by first using the ODP program to form 
the normal equations on an arc by arc basis for up to 250 parameters per arc 
selected from the following* 

« Geopotential coefficients through C 30,30 and S 30>30 

• Tracking station coordinates 

• State vector 

a Drag and solar radiation constants. 

These normal equations can be formed with or without the use of independent 
a-priori estimates and statistics. The resulting set of normal equations can 
then be combined to form a multiple arc solution for up to 500 parameters by 
using the MERGE and SOLVE programs. As m GEOSTAR-I, the EIGENVALUE 
program can be used to analyze die correlations between the elements of the 
parameter set appearing in a given system of normal equations. 

An additional program, ORBGEN has been added to the system which allows for 
independent analysis and calibration of the integration process and the orbital 
element sensitivity to particular geopotential coefficients. ORBGEN consists of 
a modified version of the ODP integration subprogram, ORBIT, together with a 
driver and subroutine which compute integration performance statistics, such as 
average stepsize and number of force model evaluations. 

The operatmg method for the orbit generator program is essentially the same 
as that used for the ODP and is documented in Appendix Al. 

III. METHOD 

3.1 Geopotential Coefficient Estimation 

Essentially, the geopotential coefficients are estimated by augmenting the exist- 
ing GEOSTAR I ODP normal equations. This augmented system of normal equa- 
tions is solved within either the framework of the single arc differential correc- 
tion process or the multiple arc processing using the MERGE/SOLVE programs. 

The normal equations in the GEOSTAR-II ODP are of the form 

[A T WA + P" 1 ] Ax = A* WA0 ~ P" 1 (2 Ax) (1) 
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where if m is the number of observations, and n the number of parameters, then 

A - an m x n matrix of measurement partials with respect to the parameters 
to be estimated, i.e., A = { a i }, where, if M l is the i th observation and 
x^ is the j th parameter to be estimated, then 

l — 1) 2, • * m 


J - 1, 2, n 

W = diagonal matrix of measurement weights 
P Q = a-priori covariance matrix on the a-priori estimates of the parameters 
Ax = n x 1 parameter correction vector to be determined 
AO = m x 1 residual vector 

2 Ax = the accumulated correction vectors over previous iterations, starting 
with the a-priori estimates. 

In the ODP program, this system can be formed with respect to the satellite state, 
up to 75 measurement station position coordinates, satellite drag and emxssivity 
factors, measurement biases, and up to 244 geopotential coefficients with a 
maximum of 250 parameters total. 

Matrix equations of the form (1) are then either directly inverted or processed 
using the MERGE and SOLVE programs for multiple arc solutions as in the 
GEOSTAR-I system See reference 1 for a detailed description of this process. 

3.2 Numerical Integration Starting Procedure 

The numerical integration of the equations of motion for position and the associ- 
ated variational equations for position partials is performed by using a Cowell 
mtegration process identical to that used in GEOSTAR-I. The Runge-Kutta | 
type starting procedure however, has been replaced by a more efficient iterative 
scheme. This process is based on closed form Newtonian-type mtegration 
formulas given by 


3M 

a = I 

13 3x 


x = h 2 
] 


n s n + o 


^ (j) x *- 1 
1 = 0 
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K o) x 



where 

x , x , x — position, velocity and acceleration vectors or vector partials with 
respect to a parametric being estimated, evaluated at time t = t n , 

h j= the stepsize to be used in the integration process, 

k ■= p - 2, where p is the order of integration process, k x j= 10, 

j = an integer which determines which time point is being computed, 

*S . ng - the first and second sums of the acceleration vectors, which, can be 

n n 

defined as 


^ = V"1 X 
n n 


n S = V" 1 *S = V" 2 x , 

n n n 

a Cj"), a* (j ) - integration coefficients which are functions of j and p. 

The coefficients and a* can be defined by integrating the appropriate New- 
tonian interpolation polynomial. The results can be expressed as 

m.) =(.»*£ (3 ^ i(J , 

'is ~ 1 


K. 

a* (]) = (- 1 y 2 ^ 



0) 


where we have the recursions 


A (3) = (- 1 Y ^ ^-i 0) 
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P (0) = - 2 rrrn % m 

-t= o 

1 

-C-o 

1 

^<o>=2>i<o^:-.t (0) 

^ = 0 


/3'/ (A 


r n «» ^ (si 

-U = 0 


and ( ) are binomial coefficients. These formulas are derived m reference 2, 
see also reference 3. 

The starting algorithm is to obtain the values 5^ , x tt , \ , i- ±1, ±2, ±3, ±4, ±5 by 
successive approximations. Letting x x ' denote the ^ th approximation, the (Z + l) st 
approximation is given by the followmg procedure: 

(i) compute the sums x S s , n S 5 usmg 



i = 0 




X. 

5 h 2 
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-L 

1 = 0 


a 




(li) Compute position and velocity vectors 


St> + O _ 


h 2 


+ 


6) x S s + 


/ °h (x - 5) x 5 - x j 


1 = 0 


5 





= h 


IS 5 + 


7^ < (1 - 5) 


i “ 0 


1 = +1, ±2, 43, ±4, ±5 


(m) Compute accelerations + 1 ^ using the force model, eompletling on 
iteration. Steps (l) - (m) are repeated until convergence. 

For the satellite position and velocity, the 1st approximation (-6=1) is obtained 
by converting the epoch position and velocity vectors to Keplenan elements and 
computing a 2-body solution at the required t x , i=±l, ±2, ±3, ±4, ±5 , For the 
variational equations the values of the position and velocity partials at t x are 
taken to be equal to their values at t 0 . 

IV. NEW AND UPDATED MODULES IN GEOSTAR-I SYSTEM 


This section describes the modification and additions made to the subroutme 
structures of the GEOSTAR-I System. 


To obtam the GEOSTAR-II ODP the GEOSTAR-I ODP was modified to allow for 
an increase m the number of geopotential coefficients to be estimated as well 
as the number of station position coordinates and measurement biases. To ac- 
complish this, the following modules were modified- 


MAIN 
OUTPUT 
OPT CRD 
STATRD 
BLKSTA 


COEFF 

DRAG 

EGRAV 

GDET 

INPUT-BLOCK DATA 


MMATRX S QUANT 

ORB 1 SWTEST 

PREDCT VEVAL 

READGP ESTIM 

SOLVGP 


Additionally, the following modules were developed specifically for the GEO- 
STAR-H ODP: 


r 


For 

Parameter 

Estimation 


< 


WRITEB 

READAP 

SYMINV 


r. 


For j 
Numerical < 
Integration 


IBINC 

MSCOEF 

MSTART 

TRANS 
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To achieve the capability of being able to integrate variational equations for a 
total of 250 parameters, the following existing GEOSTAR-I orbit generator 
modules were modified m order to efficiently utilize the space required. Those 
modules modified are: 

ORBIT HEMINT 

CSTEP SUMS 

FRCS 

In the following sections, those modules listed above which are either new or 
significantly modified existing subroutines will be documented in detail. The 
remaining modules which received min or modifications are only briefly out- 
lined, with changes indicated. Further details of these, as well as unmodified 
GEOSTAR-I ODP subroutine can be found m Reference 1. 

4.1 Modifications to Existing GEOSTAR-I ODP modules 

The modifications made to the GEOSTAR-I ODP executive routine, MAIN, and 
some of the other existing modules are designed to: 

• Call the new subroutines as listed above m proper sequence 

• Extend the program options applicable to geopotential coefficients and 
integration features 

• Extend total number of parameters to be estimated to a maximum of 
250 

• Extend tables to provide space for storing the position partials necessary 
for parameter e stim ation 

• Sort geopotential coefficients and station coordinates, being estimated 
mto an order readily acceptable to the SOLVE program 

• Extend program options to allow a-priori values of the drag and solar 
radiation constants to be different from their respective estimated 
values 

• Keep only the upper triangular portion of the normal equations matrix 
in core and invert this matrix when requested. 


Note WRITEB replaces WTMAT, TRANS replaces POSVEL, and SYMIN replaces DNVERT m the 
GEOSTAR-I ODP version 
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A summary of these modifications for the GEOSTAR-II ODP follows 


• MAIN 


DRAG "" 
EGRAY 
INPUT 1 

• MMATRX j 
SWTEST 
VEVAL _/ 

PREDCT 1 

• GDET > 
ORB 1 J 

BLKSTA 'J 

• COEFF ^ 
SQUANT J 

e STATED 


e OUTPUT 


• OPT CRD 


• READGP 


• ESTIM 


- the ODP control program. Modified to control the new 
subroutines WRITEB, READAP 


Modified to extend tables for storing the position 
- partials with respect to the geopotential coeffients 
being estimated 


Modified to accommodate the additional measurement 

- partials needed for the additional number of parameters 
being estimated 

- Modified to provide space for additional station position 
coordinates bemg estimated 

Modified to extend the total number of station position 
coordinates bemg estimated to a maximum of 75 

_ Modified to include the I/O required for the a priori 
values of the drag and solar radiation constants 

Modified to sort* the station position coordinates into 

- the proper order for easy access when creating a B 
matrix 

Modified to sort* the geopotential coefficients being 

- estimated into the proper order for easy access when 
creating a B matrix 

Modified to keep only the upper triangular portion of the 

- normal equations matrix, and to allow a full variance- 
covariance matrix for state 'only 


*Sorting was previously done in the subroutine WTBMAT of the GEOSTAR-I ODP 
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ORBIT'} 
CSTEP I 
* FRCS [ 
SUMS J 


Modified to utilize efficiently the space required for the 
- integration of the equation of motion and the variational 
equations 


• HEMINT 


Modified to interpolate for the position partials more 
efficiently. 


4.2 New GEOSTAR-II Modules and Significantly Modified GEOSTAR-I ODP 
Modules 


The following GEOSTAR-II ODP modules were written to allow for the use of a- 
priori estimates for the drag constant, solar radiation constant, and geopotential 
coefficients, I/O interface requirements with the LUNGFESH MERGE and SOLVE 
1 programs; and numerical integration of the equations of motion. 

The starting procedure portion of the orbit generator portion of the GEOSTAR-I 
ODP program was replaced by a new set of subroutines which improves the 
efficiency of the system. All other routines which comprise the orbit generator 
portion of the system were modified for more efficient core usage. 

Those subroutines effected by the new starting procedure are summarized in 
Table 1. 


Table 1 


Replaced GEOSTAR-I ODP Routmes 

New GEOSTAR-H ODP 
Routmes (approximate equivalents) 

TABLE, RK 

IBINC, MSCOEF, MSTART 
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ENTRY 0 


Purpose To retrieve the entry point address (EPA) at which the 
program is loaded mto core at execution time. 

Called By: 

MAIN 
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Calling Sequence: 


CALL ENTRY 0 <KKK1) 


Variables Not m Common: 

Fortran Name Format Description 

KKK1 I Entry point address expressed in 

decimal. This information is passed 
by COMMON/ENTRY/ to subroutine 
OUTPUT where it is printed. 
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ESTIM 


Purpose: Sums all measurement partials mto the normal equations matrix 
and right hand sides for each observation data point. When all 
observations have been processed a-priori parameter weights 
(input to the system as parameter sigmas or a full weight 
matrix) are applied. 


Called By: 


MAIN 
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Method: 


The normal equations are formed by summing all measurement partials for each 
observation data point and storing these values into a square matrix. However, 
since the normal matrix is symmetric, only the upper triangular portion of the 
matrix needs to be saved. Taking this fact into consideration, it is possible to 
reduce the matrix storage requirements from n 2 locations ton (n + l)/2 loca- 
tions (where n represents the number of parameters being estimated) by storing 
the equivalent symmetric two dimensional array into a single dimensioned array. 
Corresponding array locations mthe single dimensioned array SI,* are com- 
puted from the function 


KI (I, JD = MATD (I - 11 + J - (I - 1) 1/2 (1) 

where J > I and MATD is the total number of parameters to be estimated. 

In addition to this change, the subroutine was further modified so that a non- 
diagonal a-priori variance - covariance matrix is allowed only for the orbital 
elements.** 


*The SI single dimensioned array corresponds to the two dimensional array SUM! in the GEOSTAR-I 

ODP 

**GEOSTAR-l ODP on option allowed a full variance - covariance matrix for all the parameters 
being estimated 
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Calling Sequence 


CALLESTIM(ITYPE,SIGMA. s RESID, PMPXO, PMPSTA, ISTA1, INDXBS, 
INDXTM, TPART, MATD) 

COMMON Blacks Used. 


CONST1 

NON2 

CELEM 

SETSW 

CQUANT 

SEQ 

PRIORI 

HELP 

ESTGP 

NON4 


Variables Not m COMMON: 


FORTRAN Name 

Format 

Description 

ITYPE 

I 

1 ■= Initialization, 2 = Increment, 3 = Add a-priori 
information 

SIGMA 

R 

Measurement sigma 

RE SID 

R 

Measurement residual 

PMPXO (8) 

D 

Array of measurement partials with respect to 
satellite parameters (BM/Bx o ) 

PM PST A (3) 

D 

Array of measurement partials with respect to 
station position (B M/ BS) 

ISTA1 

I 

Total number of stations to be estimated in initial- 
ization 

INDXBS 

I 

Array index for measurement bias 

INDXTM 

I 

Array index for measurement timing bias 

TPART 

R 

Measurement partial with respect to timing bias 
(BM/BB) 

MATD 

I 

Total number of parameters to be estimated. 
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FORTRAN Name Format 


Description 


NPAR 


I Number of satellite parameters to be estimated 


NSTA 


I Total number of measurement stations to be 

estimated. 


NOEST 


True if measurement station is not to be 
estimated 
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IBINC 


Purpose: To compute the binomial coefficients needed to compute the 
coefficients for the multistep starter formulas. 

Called By: 

MSCOEF 
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Method: 


IBINC computs by addition, a (n + 1) X (n + 1) matrix, BC, containing 
Pascal's triangle. All upper triangular elements above the major diagonal 
are zero. Any bmonial coefficient up to (") can be selected from this 
matrix by the indices (i + 1 , n + 1). 
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Calling Sequence: 


CALL IBINC (N, BC) 

Variables Not m COMMON: 

FORTRAN Name Format Des cription 


BC (n + 1, n + 1) 


D 


N 


I 


Binomial 
coefficient order 



MSCOEF 


Purpose. To compute the coefficients needed for the multistep starter 
formulas. 

Called By: 

MSTART 

Calls 

IBINC 
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Method 


The coefficients a , a* for the formulas 

1 9 i- 


x ( t + s h) = h 2 


V 2 x (t) + (s - 1) V" 1 x (t) 


n 


+ 



a ± (s) x (t - ih) 


x ( t + s h) 


= h 




a* (s)x(t-ih) 


are computed, for any s, by the following: 


a 

1 


(s) = (- 1 y 



(s) 


< (s) = (- 1 y 



i = 0, 1, 2, • n 


where (^) are the bmominal coefficients and where the fi" (s), /?' (s) are ob- 
tamed by the formulas 


P[ Cs) 



-4 ^ 


/S" Cs) 


l 

2] % CO) A _e C si 

"£•= o 
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where 


jS (s) =(- 1 Y 


P[ ( 0 ) = 


i - 1 



0" (0) = 2^ (°) ^-4 CO) 

-£= o 


i = l, 2, 


and 


P 0 (0) = /?'' CO) = /S' (s) = /3'' (s) 
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Calling Sequence: 


CALL MSCOEF (M,N) 


COMMON Block Used: 

MCOEFF 

Variables Not m COMMON- 
FORTRAN Name Format 

Description 

M 


I 

M = -s 

N 


I 

Number of coefficients to be generated, 

BPZERO (20) 


D 

/?' (0) 

BIS (20) 


D 

£ Cs) 

BPS (20) 


D 

(s) 

BDPZ (20) 


D 

j8" (0) 

BDPS (20) 

( 

D 

/3" (s) 

BINCI (16, 16) 


D 

Bmomial coefficients 
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MSTART 


Purpose: To compute the required starting values for the integration of 
the equations of motion and the variational equations using an 
iterative procedure based on closed form multistep formulas. 


Called By. 


ORBIT 


Calls 


TRANS 
FRCS 
VEVAL 
MS COE F 
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Method 


An iterative process based on the closed form of the Newtonian-type integration 
formulas is used to compute the required starting arrays of accelerations and 
acceleration partials. These formulas are given by: 


n + } 


= h 2 


<i 



a 


) 

0 ) * , 


where 


’ +j ~ 



0 ) 





(1) 

( 2 ) 



n n 


Initial first approximations for , x* 1 ) , i = ±1, 2, 3, 4, 5 are obtained 

using a 2-body solution. Next, the first and second sums, iS g , i© are computed 
and used in equations <1) and (2) to obtain corrected values for position and 
velocity, x(^ +1 ), x Corrected value of the accelerations are computed 

whenever a successive correction is necessary. 
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Calling Sequence: 


CALL MSTART (IEQ) 

COMMON Blocks Used: 

WORKER 

ANPART 

GRBLOK 

XYZ 

MCOEFF 

LIMITS 

Variables Not m COMMON: 

FORTRAN Name Format Description 


M 

I 

k 

MN 

I 

k/2 

TX (3) 

D 

Components of position vector 

TXD (3) 

D 

Components of velocity vector 

TR 

D 

Time of interest 

ITER 

I 

Number of iterations needed to produce starting 
arrays 

MNI 

I 

M-I __ 

Si (3, 50) 

D 

First sum array 

S2 (3, 50) 

D 

Second sum array 

SCDX (3) 

D 

Corrected value of position or position partials 

SCOXD 

D 

Corrected value of velocity or velocity partials 

SUMI 

D 

Intermediate quantity in computation of position 


vectors 
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FORMAT Name 

Format 

Description 

SUM2 

D 

Intermediate quantity m computation of velocity 
vectors 

REAR 

D 

Square root of the sum of the squares of the 
components of the position vectors 

SUM 

D 

Maximum value of the corrector components for 
the position vectors 

SUM3 

D 

Maximum value of SUN/RBAR for all the position 
vectors 

EST (15) 

D 

Largest correction component of all equations for 
each starting position vector produced 

SXDD (3 , 20) 

D 

Saved values of acceleration vectors for all 
equations 
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REA DAP 


Purpose: To read in data to modify the prestored a-priori values and 
- stored Si g- ma va i ue s f or the geopotential coefficients 

Called By: 

OPT CRD 
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Method 


The call to READAP is initiated by the APRIOR option card. Additional cards 
which contain the coefficient name and/or an a-priori value and/or a sigma 
value are used to replace the prestored data in COMMON/ESTGP/. 
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Calling Sequence: 

CALL READAP (INTP, IARRAY, DTI, T2, T3, T4, T5) 

COMMON Blocks Used 
ESTGP 

Variables Not in COMMON. 

FORTRAN Name Format Description 

INTP I INPUT set, usually FT05 

IRRAY (4) I Not used 

DTI, T2, T3, T4, T5 D Not used 

VAL R A-priori value for the geopotential coefficient 

SIG R Sigma value for the geopotential coefficient 
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SYMINV 


Purpose: To invert a symmetric matrix stored m a single dimensional 
array. 

Called by: 

SOLYGP 


30 



Method 


The Gauss -Jordan method of condensation with partial (column) pivoting is used 
The inversion is performed on a single dimensional array containing the upper 
triangular portion of the symmetric matrix. 
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Calling Sequence- 


CALL SYM3NV (SI, NPARAM,NPA RAM, DELTA) 
Variables: 


FORTRAN Name Format 
SI (3137 5) D 


Description 

Input matrix to be inverted and output inverted 
matrix 


NPARAM 
DELTA (250) 


I Total number of unknowns to be estimated 

D Dummy array 
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TRANS 


Purpose: To convert osculating 2-body orbital elements at epoch to 
inertial position and velocity vectors at time t > to. 


Called By: 


MAIN 

MSTART 
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Method: 


Given the Keplerian orbital elements at epoch (t = to ) 
a Q = semi -major axis 
e Q = eccentricity 

Z = mean anomaly 
o 

= argument of perigee 
h Q - argument of node 

and a time t, a 2-body position and velocity vector is computed by: 

(a) computing the mean anomaly 

l + n ( t - t Q ) 

where 



GM the gravitational constant times mass of the earth, 

(b) Using this mean anomaly Z , solve Kepler T s equation 

E-esmE - Z 


by successive approximation for the eccentric anomaly E, 

(c) Compute vectors P , Q by 

P x = cos g Q cos h 0 - sin g Q sin h Q cos I Q 

P y = sm g Q cos h Q cos I Q + cos g Q sm h Q 
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p 2 = sin g 0 sm I Q 

Q x = - sm g 0 cos h Q - cos g Q sm h Q cos I Q 
Q y = cos g Q cos h Q cos I Q - sm g Q sm h Q 
Q z = cos g Q sm I 0 

(d) Compute position and velocity vectors [x , x ] by: 
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Calling Sequence: 

CALL TRANS (T, TX, TXD) 

COMMON Blocks Used- 

SCELEM 

CONST 1 


Variables Not m COMMON: 


FORTRAN Name 

Format 

Description 

T 

D 

Request Time 

TX (3) 

D 

Rectangular components of the position 
vector 

TXD (3) 

D 

Rectangular components of the velocity 
vector 

ELEM (1) 

D 

a o 

ELEM (2) 

D 

e o 

ELEM (3) 

D 

T o 

ELEM (4) 

D 


ELEM (5) 

D 

§0 

ELEM (6) 

D 


N 

D 

3 

L 

D 

Mean anomaly 

E 

D 

E 

DE 

D 

SE 

D 

D 

n 

1 - e cos E 

P(3) 

D 

P 


36 



Fortran Name 

Format 

Description 

Q (3) 

D 

Q 

RXV 

D 

I X X X 1 
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WRITEB 


Purpose; To output on a specified device the normal equations, together 
with the latest values and the labels of the parameters, in the 
SOLVE B matrix format. 


Called By: 


MAIN 



Method: 


The parameter sequence of the normal equations in the GEOSTAR-II ODP is dif- 
ferent from that required by the SOLVE B matrix format (see Reference 1). The 
subprogram WRITEB creates a new normal equation matrix m the B matrix 
format by transferring rows and columns from the original GEOSTAR-H ODP 
matrix. Blocks of data for a particular parameter type are shifted, one at a 
time, into the designated B matrix areas. In order to conserve space, the new 
matrix in the B matrix format overlays the area previously used to store the 
integration, arrays. 

WRITEB also takes each parameter type group, and assigns labels and sorts the 
parameters into the B matrix format within each particular parameter group. 

The right hand side of the normal equations, is also in a different sequence from 
that required by the B matrix format, and is shifted into the designated arrange- 
ment. Again, to conserve space, the new right hand side m the B matrix format 
overlays the area previously used to store the first sums needed m the integrator. 
After the shifting, sorting, and labeling has taken place, the normal equations 
matrix is written by rows m the B matrix format. 

The total variance, as required by the SOLVE program, is computed as- 


V 1 = (RMS) 2 * [NOB - (NPARAM -4)] 


where 


RMS = total weighted standard deviation 
NOB = total number of observations used 
NPARAM = total number of parameters being estimated. 
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Calling Sequence: 

CALL WRITES (RMSTOT, NOWTOB, NSTEST , NBIAS, BNAME, BMT) 
COMMON Blocks Used: 

SEQ ESTGP 

HELP CELEM 

NON2 PRIORI 

STANUM UMITS 

CONST 1 WORKER 


Variables Not In Common: 


FORTRAN Name 

Format 

Description 

RMSTOT 

R 

Total weighted standard deviation 

NOWTOB 

I 

NOB, number of weighted observations used in 
computing RMSTOT 

NSTEST 

I 

Number of stations whose coordinates are to be 
estimated 

NBIAS 

I 

Number of bias parameters parameters 

BNAME (3) 

R 

B matrix name (12 EBCDIC characters) 

BMT 

I 

Output data set reference number 

VI, V2 

D 

Total variance 

NPARAM 

I 

Total number of parameters to be estimated 

BMATRX (31375) 

D 

The normal equations m the SOLVE B matrix format 

BRHS (250) 

D 

The right hand side in the SOLVE B matrix format 

SBUF (250) 

D 

Dummy storage arrays 
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V. COMMON BLOCK VARIABLE DESCRIPTION 


The following sections contain the COMMON block descriptions of the new 
COMMON areas used m the GEOSTAR-II ODP. These descriptions include the 
variables contained in these areas , their type and meaning, and the subroutines 
which either define or use these variables. Also, following the section of the 
COMMON block variable descriptions, a section containing a cross reference 
table of all the COMMON areas as used in each subroutine is presented. 

5 I Modified GEOSTAR-I ODP COMMON Blocks 

This section contains the names of the COMMON Blocks, the variable names and 
the dimensions of those variables that were changed from the GEOSTAR-I ODP 
for the GEOSTAR-II ODP. 
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COMMON/SCRTCH/ 

common/ seq/ 

COMMON/C PARTL/ 
COMMON/NONl/ 
COMMON/ NON3/ 
COMMON/NON4/ 

C OMM ON/NON6 / 


COMMON/ ESTGP/ 
COMMON/ PRIORI/ 


COMMON/ LIMITS/ 
COMMON/WORKER/ 


COMMON/ PRIOR/ 
COMMON/ LIMITS/ 

C OMMON/WORKER/ 


PSIG (3,3,25) 

SUM2 (250), PARAMS (250), LABEL (250) 

PXPXO (250,6) PMPXO (250,2) 

PLHSW (25) 

TTL (250) 

BSEND (244), BSTRT (244), BYT PE (244), 
BSTANO (244), ISTEST (25) 

XYZNOM (25,6), PLHNOM (25,6), HN (25), 

SLATN (25), LATDN (25), LATMN (25), 

LONDN (25), LONMN (25), SLONN (25) 

CSA (250), CSE (250), CSSIG (250), 

NCSN (250), NCSM (250) 

STASIG (3,3,25), PLHSIG (3,3,25), 

XSTA (25), YSTA (25), ZSTA (25), 

PBIAS (244), BIAS0 (244), BIASSG (244) 

SI (3, 250), S2 (3, 250) 

X(6 ,3,250), XD (6, 3, 250), XDD (20, 3, 250) 

XXDD (3,250), VRCOV* (14000), JHIGH*, VRCSW* 

NOTES 


Written over m CSTEP 
Written over in MAIN 

*Dummy array and dummy variables used to increase 
the size of the COMMON block. This was done because 
the subroutine WRITEB uses this area when the normal 
equations matrix is put into the B matrix format. 
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5.2 New GEOSTAR-II ODP COMMON Blocks 


This section contains a description of the new COMMON areas used in the 
GEOSTAR-n ODP program. 
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/entry/ 


COMMON/ENTRY/KKK 1 


Variable 


Program Where 

Program Where 

Name 

Format D es cnption 

Defined 

Used 

KKK1 

l Entry point address 

MAIN 

OUTPUT 


where the program 
is loaded into core 
at execution time. 
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/HELP/ 


COMMON/HE LP/S1 (31375) 


Variable 


Program Where 

Program Where 

Name Format 

Description 

Defined 

Used 

SI (31375) D 

Normal equations 

ESTIM 

MAIN 


matrix 


ESTIM 


SOLVGP 
WHITE B 
ORBIT 


MSTAKT 

During the initialization process of the integrator, the storage area m 
COMMON/ HELP/ is utilized by subroutine MSTART to store temporary integra- 
tion arrays . The subroutine ESTIM then overlays this area with the normal 
equations matrix. After solution, the matrix SI (31375) contains the variance- 
covariance matrix of all parameters which are adjusted by the differential 
correction. 
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/IPERS/ 


COMMON/IPERS/ITRS, NSTEP 


Variable 

Name Format Description 


Program Where 
Defined 


Program Where 
Used 


ITRS I Total number of ORBIT 

predictor-corrector 
iterations for equa- 
tions of motion m 
CSTEP 

NSTEP I Total number of ORBIT 

Cowell steps for 
the equations of 
motion 


CSTEP 

MAIN 


CSTEP 

MAIN 
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/MCOEFF/ 


COMMON/MCOEFF/SBPS(16), SBDPS (16), STOL, II, 12, N00 


Variable 
Name Format 

Description 

Program Where 
Defined 

SB PS (16) 

D 

Coefficients for multi- 
step starter equations to 
compute positions 

MSCOEF 

SBDPS (16) 

D 

Coefficients for multistep 
starter equations to 
compute velocities 

MSCOEF 

STOL 

D 

Convergence criteria 

BLOCK DATA 
OPT CRD 
(initialized to 
1.0 x 10" 13 ) 

11 

I 

Value zero 

MSCOEF 

12 

I 

Value zero 

MSCOEF 

N00 

I 

Total number of iter- 
ations allowed 

BLOCK DATA 
OPT CRD 
(initialized to 
15) 


Program Where 
Used 


MSTART 


MSTART 


MSTART 

OUTPUT 

MSTART 

MSTART 

MSTART 
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/SCLEM/ 


COMMON/SCELEM/ ORBEPC (6) 

Variable Program Where Program Where 

Name Format Description Defined Used 


Orbital Elements - Osculating 
Keplenan 


ORBEPC (1) 

D 

Semi -major axis- 
meters 

MAIN 

TRANS 

ORBEPC (2) 

D 

Eccentricity 

MAIN 

TRANS 

ORBEPC (3) 

D 

Inclmation-radians 

MAIN 

TRANS 

ORBEPC (4) 

D 

Longitude of ascending 
node-radians 

MAIN 

TRANS 

ORBEPC <5) 

D 

Argument of percenter- 
radians 

MAIN 

TRANS 

ORBEPC (6) 

D 

Mean anomaly-radians 

MAIN 

TRANS 
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/SEQ/* 


C0MM0N/SEQ/SUM2 (250) , PARAMS (250), IA BE 1.(250) , IDMAT 


Variable 

Type 

Description 

Program Where 
Defined 

Program Wh 
Used 

SUM2(250) 

D 

Right hand side of normal 
equations in GEOSTAR-II 
ODP parameter order 

ESTIM 

ESTIM 
WRITE B 
SOLVGP 

PARAMS 

D 

Parameter values m B 
matrix order by row 

WRITE B 

WRITE B 

LABEL(250) 

I 

Numberical identifiers for 
the type of parameters m 
PARAMS, in the same 
order 

WRITE B 

READGP 

WRITEB 

IDMAT 

I 

B matrix identification 
number. 

OPTCRD 

WRITEB 

OUTPUT 


* Replaces COMMON/BE Q/of the GEOSTAR- 1 ODP 



/SINIT/ 


COMMON/SINIT/ITERSW 


Variable 

Type 

Description 

Program where 
Defined 

Program where 
Used 

ITERSW 

L 

Switch to ini- 
tialize the 
starter for 
each iteration 

ORBIT 

MSTART 
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5.3 GEOSTAR-H ODP COMMON Block Cross Reference Table 


This section contams a cross reference table describing the COMMON area 
structure m the GEOSTAR-II ODP. 
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COMMON BLOCKS 



COMMON BLOCKS 


GEOSTAR-I COMMON Block Cross Reference Table (continued) 


SUBROUTINES 


PPPPRRRRSSSS 

LRRREEEEOQTT 
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S 

L 
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X 

R 

X 

R 
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V 

T 



B 
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T 


ABCOEF 














X 










ANPART 


















X 






CELEM 



















X 





CGEOS 








X 




X 












COFIT 


X 
















X 






CONST1 

X 

X 








X 


X 





X 

X 

X 





CONST2 


X 










X 












CONST3 


X 


X 











X 



X 






CONVRG 












X 












CPARTL 


X 






















CQUANT 

X 

X 





X 



X 


X 












CSTHET 


X 


















X 

X 

X 

X 

CSTINF 











X 













CSVEVL 













X 





X 






DRGBLK 


















X 






ESTGP 


X 



X 

X 






X 




X 


X 

X 





FLXBLK 


X 
















X 






FMODEL 






X 







X 











GRBLOK 


















X 






HELP 









X 










X 





LIMITS 


X 










X 




X 


X 

X 





MOONGR 















X 



X 






NON1 












X 












NON2 


X 




X 






X 







X 
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X 
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- 








X 
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X 












NON6 












X 
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X 
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X 




















PRIORI 









X 

X 


X 







X 





SCELEM 

















X 







SCRTCH 










X 


X 












SEQ 






X 



X 










X 





SETSW 
















X 


X 






S1GMAC 












X 












STANUM 

X 






X 





X 







X 





VRBLOK 


















X 






WORKER 














X 


X 

X 


X 





XY2 


















X 






XYZOUT 


X 


X 
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VI. PROGRAM FLOW CHARTS AND SUBROUTINE SUMMARIES 

The following sections contain the flow diagrams of the mam control routmes 
in the GEOSTAR-IE ODP and summary of new modules used in the program. 

Cross reference tables detailing the entire subroutine structure of the GEOSTAR- 
II ODP is also presented. 
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6.1 GEOSTAR-II QDP Flow Diagrams 

This section contains the flow diagrams of the GEOSTAR-II ODP executive 
program MAIN and the orbit generator control subprogram ORBIT. 
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MAIN PROGRAM 
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ORBIT SUBROUTINE 


C 


ORBIT 


INPUT 
REQUEST 
TIME t D 


3 



INITIALIZE 
CONSTANTS 
FOR ORBIT 
GENERATOR 
t= t„ 


EPHQAN 


COMPUTE 

EPHEMERIS 

QUANTITIES 


T 


MSTART 


COMPUTE 

STARTING 

VALUES 


T 


SUMS 


COMPUTE 

INITIAL 

1st AND 2nd SUMS 


I 



NO 


1 

YES 

UPDATE 

INTEGRATION 

TIME t 

= t + h 

’ 



CSTEP 


INTEGRATE THE 
EQUATIONS OF MOTION 
AND REQUIRED 
, VARIATIONAL EQUATIONS 
TO TIME t 


AA 



HEM1NT 


INTERPOLATE FOR 
POSITION, VELOCITY ANDJ 
PARTI ALS AT TIME 


c 


RETURN 


3 


AA 
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6.2 GEOSTAR-n ODP Subroutine Summary 

This section contains a brief sum m ary of the new modules used in the GEOSTAR- 
II ODP. Complete documentation of the GEOSTAR -I ODP modules can be found 
m Reference 1. 
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Subroutine Summary 


ENTRYO 

IBINC 

MSCOEE 

MSTART 

READAP 

SYMINV 

TRANS 


Retrieves the entry point address (EPA) at which the program is 
loaded into core at execution time. 

Computes the binomial coefficients needed to compute the co- 
efficient for the multistep starter formulas. 

Computes the coefficients needed for the multistep starter 
formulas . 

Computes the required starting values for the integration of the 
equations of motion and the variational equations using an iter- 
ative procedure based on closed form multistep formulas. 

Modifies prestored a-priori and sigma values for the geopotential 
coefficients. 

Inverts a symmetric matrix which is stored m a smgle dimen- 
sioned array. No restrictions on the dimension of the input 
matrix. 

Converts osculating 2~body orbital elements at epoch to inertial 
position and velocity vector at time, t > t Q . 
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6.3 GEOSTAR- II ODP Subroutine Cross Reference Tables 

This section contains cross reference tables detailing the subroutine structure 
m the GEOSTAR-II ODP executive program MAIN and the orbit generator con- 
trol subprogram ORBIT. 
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/ 

MAIN Program Cross Reference Table 


CALLING SUBROUTINES 

MSCC ODD DDDD E E E F O M M M O O O 

A LOS aae IJOR LPQRESSUBPR 
1 KETTYN FUDA EHUCOCTTSTB 

ns F5 EE§ Fb§c M Q A S S O A A D C 1 

TFPSAT n AT REkTOR 

A RY D NR DFTETD 


w 


a 


a 


u 


o 


ADDYMD 

AffiK 

APPER 

COEFF 

COEFL 

DARCTN 

DATES 

DA YEAR 

DEKWL 

DENSITY 

DIFF 

DIN RAD 

DJUL 

DODSRD 

DOT PRD 

DRAG 

EGRAV 

ELEM 

ELEMK 

ENTRY© 

EQN 

EQUATR 

ERROR 

ESTIM 

FRCS 

GDET 

Ct'OSRD 

GTIMtN 

GTIMOT 

HEM I NT 

101 NC 
1NV2 
1NV3 

MATRIX 

MOO NAD 

MSCOEt 

MULMAT 

NUMDFR 

NUTATE 

OBSDOT 

OPTCRO 

ORBI 

ORBIT 

OUTPUT 

PRECES 

PREDCT 

HLADAP 

READGP 

REFCOR 

ROTMAT 

1IYMDI 

SATCL2 

SATCU. 

SDOLL- 

SDTPRD 

SOLVCP 

SQU ANT 

STAINF 

STATEUI 

STORGP 

SUN 

SUNCRV 

SYMMET 

TD1F 

TRANS 

VEVAL 

WR1TEB 

XEF1X 

XINERT 

YEFIX 

YINERT 

YMDAY 
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CALLED SUBROUTINES 
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VII. ORBGEN 


7.1 Introduction 


This section describes the modifications and additions made to the subroutine 
structures of the variable stepsize venson of GEOSTAR -I ODP m order to 
create a stand-alone orbit generator program, ORBGEN. 

The ORBGEN program can optionally numerically integrate the variational 
equations as well as the equations of motion. The limits on the types of vari- 
ational equations are as follows: 

(a) all parameters - 20 

(b) gravity coefficients - 20 

(c) state parameters - 6 

(d) drag parameters - 1 

(e) solar radiation parameters - 1. 

The partials are available with respect to either rectangular or Keplerian ele- 
ments. Also, osculating orbital elements are available at each output step. 

The program can operate m either a fixed or variable stepsize mode. In the vary 
step mode, a resume of stepsize distribution over the integrated trajectory will 
be output. In the fixed step mode, different stepsizes can be used for the equa- 
tions of motion and variational equations. Different orders can be used m either 
mode. 

Various force model options are available for the integrations. Using the elliptic 
motion option, accumulated error estimates are computed by comparing with 
an analytic 2-body trajectory. 

7.2 New and Updated Modules m GEOSTAR-I ODP 

To obtain ORBGEN the following GEOSTAR-I ODP modules were modified* 

MAIN CSTEP VEVAL 

OPTCRD TEST MSTART 

ORBIT FRCS 

In addition, the following modules were developed specifically for ORBGEN: 

RESUME 

WRITER 

APART 
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In the following sections, those modules listed above which are new, or sig- 
nificantly modified existing subroutines well be documented in detail. The re- 
maining modules which received only minor modifications are briefly outlined, 
with changes indicated. 

7.3 Modifications to Existing GEQSTAR-I ODP Modules 


The modifications made to the GEOSTAR-I ODP are designed to 


© Call the new subroutines 

© Retain only the subroutines, logic, and output pertinent to an orbit 
generator program 

• Integrate m a UTC time system 


A summary of these modules, as modified for ORGBEN follows: 


© MAIN The program control module. Modified to control only the 

integration of an orbit and associated variational equations. 


• OPTCRD Modified to mclude an option to integrate in the A1 time 
system 


a 


ORBIT 

CSTEP 

TEST 

FRCS 

VEVAL 

MSTART 


> 


j 


Modified to save various quantities which are used to compute 
a statistical summary at the end of a run 


7.4 New ORBGEN Modules 


These modules were written to allow for the computation of partials with respect 
to Keplerian elements as well as to accumulate a statistical summary of the 
integration performance 
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APART 


Purpose: To compute the Jacobian of the osculating Keplerian elements 
with respect to rectangular coordinates. 

Called by: 

WRITER 

Calls: 

INV2 
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Method 


The Jacobian. 


3 (a, e, a, fi, co, M) 
3 (x, y, z, x, y, z) 

is computed by inverting the 6X6 matrix 


3 (x, y, z) 
3 (a, e, a) 

3 (x, y, z) 
3 (a, e, i) 


3 (x, y, z~) 
3 (Q, oo, M) 

3 (x, y, z) 
3 (U <z>, M) 


The required partials can be derived from the relations 


x = R xq (Q, co, co) q (a, e, M) 
x = R xq (fi, i, co) q (a, e, M) 


where 


a (cos E - e) 

a a 

3 q = 

- sin E 

a Vl - e 2 sin E, 

Vl - e 2 cos E 

1 - e cos E 

0 


0 


and R xq is the rotation matrix which transforms the [x, x] coordinates to [q, q] 
coordinates. Explicit expressions for R and the partials with respect to the 
elements can be found in W. Kaula, "Theory of Satellite Geodesy,’' 1966, pp. 
67-68. 
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COMMON Blocks Used: 


DC 

XYZOUT 

CONST1 


Variables not m COMMON* 


FORTRAN Name 

Format 

CSI 

D 

SSI 

D 

SA 

D 

SN 

D 

OME2 

D 

SE 

D 

SCOMEG 

D 

CCOMEG 

D 

SOMEG 

D 

COMEG 

D 

B(I, J, 1) 

D 


Description 


Cos I 
Sin I 
a 

n, mean motion 
(1 - e*) 
e 

Smfl 
Cos Q 
Sin co 
Cos co 

9 (x, y, z, x, y, z ) 
B (a, e, 1 , fl, co, M) 
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RESUME 


Purpose. To write a summary of the integration statistics. 
Called by 


MAIN 

ORBIT 
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Method: 


The routine computes the integration statistics by accumulating the desired 
variables withm the subprograms where these variables are being determined 
and passes the accumulated values thru COMMON/TIMER/. 

The variable step statistics are computed by allowing a possibility of 10 step- 
size ranges of (0 - A, A - .8, .8 - 1.6 mm., etc.) and accumulating in each 
particular range the amount of orbit time spent, and the number of stepsizes 
taken within that range. The average stepsize is then computed by adding all 
the stepsizes of a particular range and dividing by the total number of step- 
sizes m that range. The percent of orbit time is computed by adding the orbit 
time accumulated and dividing by the total run time. The time that the program 
spends in setting up the initial starting table is not added into the stepsize 
analysis statistics. 

The routine has 3 entries ; 

(1) initialize phase for accumulation purposes 

(2) accumulation of stepsize statistics 

(3) final computation and output of statistical summary. 
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Calling Sequence: 


CALL RESUME (END) 
COMMON Blocks Used: 


TIMER WORKER 

LIMITS COFIT 

SSTEP IORBIT 

COWS 

Variables not in COMMON: 


FORTRAN Name 

Format 

Description 

IND 

I 

Entry point indica- 
tor 

STEPN 

R 

Old stepsize when 
a stepsize change 
occurs (sec.) 

STMES 

R 

Last time at which 
a stepsize change 
occurred (sec.) 

IRANGE(ll) 

I 

Number of step- 
sizes within a par- 
ticular range 

STEP (11) 

R 

Accumulation of 
stepsizes within a 
particular range 
(min.) 

TIME 1 (11) 

R 

Accumulation of 
orbit time within a 
particular range 
(sec.) 


STEP 


R 
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Old steps lze when 
a stepsize change 
occurs (min.) 



FORTRAN Name 

Format 

Description 

ST1 

R 

Orbit time spent at 
a particular step- 
size range 

AVITER 

R 

Average number of 
predictor-corrector 
iterations during a 
run 

DAYEND 

D 

Total run time m 
seconds 

TIMES 

R 

Percentage of the 
total orbit time 
spent in a particular 
steps lze range. 
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WRITER 


Purpose. To write the values of the satellite position and velocity vectors 
and optionally do the following: 

(1) write the osculating Keplerian elements 

(2) write the partials of position and velocity with respect to 
parameters being estimated 

(3) compute and write the partials of the osculating Keplerian 
elements with respect to parameters being estimated. 


Called By: 


MAIN 


Calls: 


TRANS 

ELEMK 

APART 
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Method: 


If the drag model is being used in the acceleration computations, the norm of the 
drag vector is automatically computed The norm is taken to be the largest 
component of the drag vector. 

If the 2-body motion option is being used, the norm of the error vectors is 
computed automatically. The analytic solution at the requested time is com- 
puted and subtracted from the position and velocity vectors computed from the 
integration; the norm is taken to be the largest component of this vector. 

The partials of the osculating Keplerian elements with respect to parameters, 
if requested, is computed by multiplying the Jacobian 

3 (a, e, 1 , fi, co, M ) 

3 (x, y, z, x, y, z) 

by the matrix 

3 (x, y, z, x, y, z) 

^ (x Q , y Q , z Q , x 0 , y Q , z Q , C nm , S nm , * * • ) 
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Calling Sequence: 


CALL WRITER (PXPXO) 
COMMON Blocks Used* 


WORKER 

XYZ 

XYZ OUT 


LIMITS 

NON3 

CONST2 


DRGBLK 

IORBIT 

COFIT 

COWS 

TIMER 

DC 


Variables not in COMMON: 

FORTRAN Name 

Format 


Description 

TREQ 

D 


Request time (time 
from epoch m 
seconds) 

CON 

D 


Cross product mag- 
nitude of position 
and velocity vectors 

DROG(3) 

D 


Component of drag 
vector 

DN 

D 


Norm of drag vector 

TI 

D 


Request time (time 
from epoch, in secs.) 

TX(3) 

D 


Analytic 2-body 
value of position 
vector 

TXD(3) 

D 


Analytic 2-body 
value of velocity 
vector 

EV 

D 


Norm of error vec- ■ 
tor for position 

EVD 

D 


Norm of error vec- 
tor for velocity 
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7.5 ORBGEN COMMON Block Variable Description 


The following section contains the COMMON block descriptions of the new 
COMMON areas used m the ORGBEN program. 
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/OPT/ 


C OMMON/OPT/OSTE P 

Variable Name Format 
OSTEP D 


Program Where 
Description Defined 


Optimum stepsize TEST 


Program Where 
Used 


MSTART 
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/TIMER/ 


COMMON/T3MER/TOTIME, TABLET(2), CTIME(2), STIMED, STEPD, NOFT{2), 
NOFC(2), N0C0WL(2), ITERS, ITERMS(2)» N0STEP(2) 


Variable 


Program Where 

Program Where 

Name 

Format 

Description 

Defined 

Used 

TOTIME 

R 

Total run time 

MAIN 

ORBIT OSTEP 

TABLET 

R 

TABLET(l) - Time spent 

MSTART 

RESUME 

(2) 


xn TABLE for equations of 
motion 

TAB LET (2) - Time spent 
in TABLE for variational 
equations 



CTIME 

R 

CTIME (1) - Time spent in 

CSTEP 

RESUME 

(2) 


CSTEP for equations of 
motion 

CTIME (2) - Time spent 
m CSTEP for variational 
equations 



STIMED 

R 

Time spent at a particu- 
lar stepsize 

TEST 

RESUME 

STEPD 

R 

Old stepsize 

TEST 

RESUME 

NOET(2) 

1 

NOFT(l) - Number of times 
FRCS routine called by 

FRCS 

RESUME 



equations of motion when 
computing starting table 
NOFT(2) - Number of times 
VEVAL routme called by 
variational equations when 
computing starting table 

VEVAL 


NOFC(2) 

I 

NOFC (1) - Number of times 

FRCS 

RESUME 



FRCS routme called by 

VEVAL 



CSTEP 

NOFC (2) - Number of times 
VEVAL routine called by 
CSTEP 
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/TIMER/(contmued) 


Variable 

Name 


NOCOWL 

(2) 


ITERS 


ITERMS 

( 2 ) 


NOSTEP 

( 2 ) 


Program Where Program Where 
Format Description Defined Used 


I NOWCOWL(l) - Number CSTEP RESUME 

of Cowell steps for 
equations of motion 
NOCOWL(2) - Number of 
Cowell steps for variational 
equations 

I Total number of predictor- 1 CSTEP RESUME 

corrector iterations for 
equations of motion in 
CSTEP 


I ITERMS (1) - Total number MSTART RESUME 

of corrector iterations 
for equations of motion 
when computing table of 
ITERMS (2) - Total number 
of corrector iterations 
for variational equations 
when computing starting 
table 

I NOSTEP(l) - Total number TEST RESUME 

of step changes when com- WRITER 

putmg starting table 
NOSTEP (2) - Total number 
of step changes when in. 

CSTEP 
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7 .6 ORBGEN Flow Charts and Subroutine Summaires 


The following section contains the flow diagrams of the MAIN or control pro- 
grams in ORBGEN, as well as a summary of new modules used m the programs. 
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MAIN: 


Controls the integration of an orbit for a particular time span to generate a satel- 
lite ephemeris at equally spaced intervals. In addition, optionally, position 
partials can be generated along with the satellite ephemeris. 
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MAIN PROGRAM (Logic Flow) 

MAiN PROGRAM 
(Logic Flow) 
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ORBIT SUBROUTINE 




INTEGRATE THE 
EQUATIONS OF MOTION 
AND REQUIRED 
VARIATIONAL EQUATIONS 
. TO TIME i 














Subroutine Summary 


APART To compute the Jacobian of the osculating Keplerian elements 

with respect to rectangular coordinates 

RESUME Writes a summary of the integration statistics for those equations 
being integrated. Some of those statistics are: average number of 
predictor-corrector iterations, total number of Cowell steps, 
number of stepsizes selected within a particular steps lze range 
and the total amount of orbit time spent in a particular steps lze 
range. 


WRITER Writes the values of the satellite position and velocity vectors 
and optionally the partials of position and velocity with respect 
, to selected parameters at the requested time. In addition, if 

osculating Keplerian elements are requested, then the osculating 
Keplerian element partials with respect to selected parameters 
will be written at the required time. 
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APPENDIX A1 


ORBGEN INPUT CARD FORMAT DESCRIPTION 

Data Cards 

The data cards used to define a case are divided into two logical groups: 

Group 1 - Cards which must be present to define each run. There are 7 of 
these cards which must appear in a fixed order at the beginning of 
the data deck. They provide labels and the orbital elements. 

Group 2 - Option cards; these cards provide options to the user such as desired 
print internal, selection of outputs and position partials. The ORBIT 
and OUTPUT cards are necessary to define an orbit generator run, The 
Group 2 cards are placed after the Group 1 cards and may appear 
in any order. 

A detailed description of the card formats for each group follows. Table 2 

defines the format .specification codes used for individual cards. The last card 

in any deck must be a standard card ending a data set (/*m columns 1 and 2). 
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Table 2. Format Specification Interpretation 

Key: n, m = integer numbers 
b -a. blank space 

Format Code Interpretation 

In digits with no decimal point right adjusted in a field of n 

col umns 

example: 25 in an 13 format: b25 

Fn, m digits with a decimal point anywhere m a field of n columns 

or digits punched with no decimal point, in which case the 
pomt will be assumed between the m-lth and mth column of 
the field 

example: 30.1 in an F5.2 format: b30.1 or 30.1b or b301b 

Dn, m digits with a decimal pomt anywhere in a field of n columns 

or digits right justified m a field of n columns with an exponent 
of the form D±XX, where XX is the power of 10 to which to 
raise the number. If there is no decimal pomt it will be as- 
sumed to be m places to the left of the "D". 
example: 30.1 m a D8.5 format: b30.1bbb or 30.1D+00 or 
301.D-01, etc. 

An n alphanumeric characters 

En, m digits with a decimal pomt anywhere in a field of n columns 

or digits right justified in a field of no columns with an ex- 
ponent of the form E±XX, where XX is the power of 10 to 
which to raise the number. If there is no decimal pomt it will 
be assumed to be m places to the left of the "E". 
example: 30.1 m a E8.5 format: b30.1bbb or 30.1E+00 or 
3 01. E -01, etc. 
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Columns 


1-72 


TITLE INFORMATION 
CARDS 1, 2, 3 

Format Description 


12A6 Identification to be printed on all output generated by 

the case 
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Columns 

1-6 


REFERENCE DATE 
CARD 4 

Format Description 


16 Reference date of the epoch elements year, month, 

day (YYMMDD) 
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EPOCH DATE 
CARD 5 


Columns 

Format 

Description 

1-6 

16 

Year, Month, Day (YYMMDD) - Epoch date* 

7-10 

14 

Hour, minute (HHMM) 

11-17 

D7.4 

Seconds (SS.SSSS) 

24-30 

13 

Satellite ID, 5-7 digit identification number. 

V 

73-80 

A8 

Satellite name to be printed m this run. 


j 

*NOTE* The epoch tune must be before the start time. 

EPOCH = date of initial position and velocity vector. 


89 



POSITION VECTOR 
CARD 6 


Columns Format Description 


Tlie X component 
The Y component 
The Z component 

or 


1-24 

D24.16 

25-48 

D24.16 

49-72 

D24.16 

1-24 

D24.16 

25-48 

D24.16 

49-72 

D24.16 

Note: 

If the absolute 


Initial position vector 
in meters 


Semi-major axis (meters)] 

Eccentricity > Keplerian elements 

Inclination (degrees) j 


Note: If the absolute value of the 2nd element on this card is less than 1.0, 


the elements on cards 3A and 3B are assumed to be Keplerian elements. 
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VELOCITY VECTOR 
CARD 7 


Columns 

Format 

Description 



1-24 

25-48 

49-72 

D24.16 

D24.16 

D24.16 

The X component! , 

, r ,1 Initial velocity vector 

The Y component > . , 

_ , ( meters /second 

The Z component 1 

or 

1-24 

D24.16 

Right ascension of ascending '' 
node (degrees) 

v 

Keplenan 

25-48 

D24.16 

Argument of perigee (degrees) 

r 

elements 

49-72 

D24.16 

Mean anomaly (degrees) 
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Group 2 - Option Cards 


All of the Group 2 heading cards have the same format, although subdivisions of 
these headings will have their own formats. Each heading card is identified by 
a name starting in column 1, which is read under an A6 format. The name is 
followed by 9 fields.* 

Field Columns Field Format 


1 

7 

11 

2 

8 

11 

3 

9 

11 

4 

10 

11 

5 

11-25 

D15.8 

6 

26-40 

D15.8 

7 

41-55 

D15.8 

8 

56-70 

D15.8 

9 

71-80 

D10.8 


*Note* Fields 5 thru 9 of the heading cards are read under a decimal format 
and MUST CONTAIN A DECIMAL POINT 

A detailed description of each option card, in alphabetical order, follows. 
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COEFG CARD GROUP 
Header Card 


Name 
COEFG P 


Columns Format Description 


1-6 A6 


7 II 


26-40 D15.8 


This card signifies that geopotential 
cards are to follow m the input stream. 
The values of the geopotential coefficients 
can be changed (prestored values are m 
BLOCK DATA) and/or position partials 
can be computed with respect to the geo- 
potential coefficients requested. 

6 = compute position partials for state 

0 = do not compute position partials for 
state 

When not zero, partials are computed 
and printed at each ephemens point.* 


*Position partials are printed unless Keplerian elements are requested at each 
ephemens point (see OUTPUT option card), in this case partials of the oscula- 
ting Keplerian are printed. 



COEFGP CARD GROUP 
EST and CHG Cards 


Name Columns 

Format 

Description 

EST 1-3 

A3 

Indicates a required position partial 
with respect to the parameter on the 
card.* 

4 

11 

0 = use the stored value 

1 = use the value on the card 

CHG 1-3 

A3 

Indicates a change m the prestored value 
for the parameter listed on the card, hut 
no position partial will be computed 

6 

A1 

Identifies this entry as a C or S 
parameter 

*Note: Either EST or CHG options 
utilized with an EST entry. 

may be used, but column 4 will only be 

7-8 

12 

Subscript N for parameter 

9-10 

12 

Subscript M for parameter 

11-25 

E15.8 

*Value for the parameter 


*Note- E formats must be right adjusted in last 4 columns of the field. 
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COEFGP Card Group 
END Card 


Name 

Columns 

Format 

Description 

END 

1-3 

A3 

Group end card which must belpresent 
when a COEFGP card is used, and 
placed directly behind array EST or 
CHG cards. 
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CTO LI 


Name Columns Format Description 


CTOLl 1-5 A6 Redefines the predictor-corrector 

convergence criterion of the integra- 
tion, CTOL 

11-25 D15.8 CTOL value; default value is 0.5D-05 

meters 

Normally the user would not use this card and thus would permit the default value 
to hold. If an error message occurs indicating failure to converge in the pre- 
dictor-corrector cycle of the integrator, a decrease m stepsize using the STEP1 
card is suggested as a corrective measure rather than a change of the CTOL 
value. 
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DRAG 


Name 

Columns 

Format 

DRAG 

1-4 

A6 


11-25 

D15.8 


26-40 

D15.8 


Description 


This card redefines the 
atmospheric drag coefficient. 

C D - satellite drag coefficient. 

Zero indication no position partials 
with respect to C D required. 

Drag perturbations will not be applied 
li C d = €. 


Note: See Appendix A2 for program determination of drag. 
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EARTH 


Name Columns Format Description 


1-5 

A6 

This card redefines the Earth’s refer- 
ence spheroid. 



The following 2 fields provide redefini- 
tion of the limits for the geopotential mod- 
el to be used, here N and M refer to the 
subscripts C n>ra and S nim . 

7-8 

211 

N+ 1, highest order to use (initialized 
to 16) 

9-10 

211 

M+ 1, highest degree to use (initialized 
to 16) 

11-25 

D15.8 

GM - Universal gravitational constant 
times mass of earth (meter /second ) 

25-40 

D15.8 

A e - Semi-major axis of Earth's refer- 
ence speroid (meters) 

41-55 

D15.8 

f - Inverse of flattening of Earth's ref- 
erence (meters) 



Default values GM = 3.9560032 D + 14 

A e = 6378165. 

F = 298 252 



These constants are available m COMMON/ 
CONST/as the variables GM1, AE1, and 
FIN VI. 

56-70 

D15.8 

- I., the standard geopotential model is 
replaced by a 2-body model m the inte- 
gration of the equations of motion 
(state) and the integration of state par- 
tials. Other forces may be included. 

71-80 

DIO. 8 

= 1., only a 2-body model is used m the 


integration of the equations of motion 
for state and the integration for state par- 
tials. No other forces are included. 


*These options override any force model specifications m the partials integra- 
tion defined by the PARTS card. 


98 



MODE 


Name 

Columns 

Format 

Description 

MODE 

1-4 

A6 

This card is used to override the 
program determined used of the vari- 
able step or fixed step mode. 


V 

11 

- 2, use variable step 

= 4, use fixed step 

This sets the internal variable, 

MODE1. 


Note: See Appendix A2 for program determination of mode and stepsize. 
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MOON 


Name 

Columns 

Format 

Description 

MOON 

1-4 

A6 

This card redefines the ratio of the mass 
of the Moon to the mass of the Earth 


11-25 

D15.8 

Value of ratio of mass of Moon to mass 
of Earth. If value is zero lunar gravita- 


tional perturbations will not be applied 

Initially the mass ratio is defined to be ' 
0.0122999 in BLOCK DATA. This information 
is available m block COMMON/ CONST/ 
as the variable MMOONI. 
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ORBIT 


Name 

Columns 

Format 

Description 


ORBIT 

1-5 

A6 

This card defines the time span of the 
orbit requested. This card must be 
present m the data deck. 


11-25 

D15.8 

Year, month, day (YYMMDD.) ~) 

Hour, minute second (HHMMSS.) j 

start 
^ print 
time 


41-25 

D15.8 

Yera, month, day (YYMMDD.) "1 

end 


56-70 

D15.8 

Hour, minute, second (HHMMSS.) ^ 

r print 
' time 


Note The time interval at which to output the satellite ephemeris is specified 
on the OUTPUT card 
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ORDER1 


Name 

ORDER1 


Columns 

Format 

Description 

1-6 

A6 

Redefines the order, ORDER(l), for the 
integration of the equations of motion 

A 

7-8 

12 

ORDER(l), default value is 11 
Maximum = 12 
Minimum = 4 
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ORDER2 


Name Columns Format Description 

0RDER2 1-6 A6 Redefines the order, ORDER(2), for the 

integration of the variational equations 
for partials of state, drag, solar radia- 
tion and the geopotential coefficients 

7-8 12 ORDER(2), default value is 7 

Maximum = 12 
Minimum = 4 
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OUTPUT 


Name 

OUTPUT 


Columns Format Description 


1-6 A6 


7 II 


10 II 


11-25 D15.8 


This card specifies the output time 
interval. This card must be present 
m data deck. 

= 1 indicates that rectangular coordinate 
ground track is requested 

= 2 indicates that ground track, Kepler- 
lan elements and/ or partials of the 
osculating Keplerian elements are re- 
quested (see COEFGP option card) 

= 3 indicates that output times are to 
be converted from UTC to A1 time. 

= 1 indicates list all non-zero gravity 
coefficients used m the geopotential 
model 

A value greater than zero m this field 
specifies that the ground track is 
referenced to a fixed time frame defined 
by the REFERENCE card. 
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PARTS 


Name 

Columns 

Format 

Description 

PARTS 

1-5 

A6 

This card redefines the force model to 
be used for the integration of the vari- 
ational .equations. It does not effect the 
force model used for integrating the 
equations of motion. 




The following options denote exceptions 
to the equations of motion force model 
for integration of the variational equations. 


11-25 

D15.8 

- 1., Exclude the standard geopot ential 
model 


26-40 

D15.8 

= 1., Exclude drag perturbations 


41-55 

D15.8 

= 1., Exclude solar radiation preturbations 


56-70 

D15.8 

= 1,, Exclude Sun perturbations 


71-80 

D15 8 

= l. y Exclude Moon perturbations 
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SAT 


Name 

Columns 

Format 

SAT 

1-3 

A6 


11-25 

D15.8 


26-40 

D15.8 


Description 

This card redefines satellite area and 
mass. 

Cross sectional area of satellite 
(meter 2 ) 

Mass of satellite (kilograms) 


Initially the satellite area and mass 
are define as 0 m BLOCK DATA. 
TMs information is available in 
block COMMON/CONST/ as the 
variables AS ATI and MSATI. 

Drag and solar radiation pressure 
perturbations will not be applied if 
either the satellite area or mass is 
equal to zero. 
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SOLRAD 


Name 

Columns 

Format\ 

Description 

SOLRAD 

1-6 

A6 

This card redefines the satellite reflect- 
ivity 


11-25 

D15.8 

S = Satellite reflectivity (default = 1.5) 

R 


26-40 

D15.8 

Zero indicates no position partials with 
respect to are required. 


Solar radiation pressure perturbations 
will not be applied if S =0. 

R 
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STEP1 


Name 

Columns 

Format 

Description 

STEP1 

1-5 

A6 

Redefines the stepsize value, CDEL(l), 

to be used for integrating the equations of motif 


11-25 

D15.8 

CDEL(l) = stepsize (secs). Default 
value is 100. seconds. 
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STEP2 


Name 

Columns 

Format 

Description 

STEP2 

1-5 

A6 

Redefines the stepsize value, CDEL(2), 
to be used for integrating the variational 
equations. 


11-25 

D15.8 

CDEL(2) = stepsize (secs.), 
Default value is 100. seconds. 



stou 


Name 

Columns 

Format 

Description 

STOLI 

1-4 

A6 

Redefines the predict or- c orr e etor con- 
vergence criterion of -the multistep 
startmg process. 


11-25 

D15.8 

STOLI value. (default .1.0 x io -13 ) 


26-40 

D15.8 

Total number of corrector iterations 
allowed in the multistep starter 
(default «= 15). 


no 



SUN 


Name 

Columns 

Format 

Description 

SUN 

1-3 

A6 

This card redefines the ratio of the mass 
of the Sun to the mass of the Earth. 


11-25 

D15.8 

Value of ratio of mass ratio If value is 


zero, solar gravitational perturbations' 
will not be applied. 


Imtally the mass ratio is defined as 
332951.3 in BLOCK DATA. This in- 
formation is available in block COMMON/ 
CONST/as the variable MSUNI. 
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TOLS 


Name 

Columns 

Format 

Description 

TOL 

1-3 

A6 

Redefines the upper and lower bounds on 
truncation error for step size modifica- 
tion 


11-25 

D15.8 

TOL1 or T , the upper bound, default 
value is .25 D-4 meters 


26-40 

D15.8 

TOL2 or T 2 , the lower bound, default 
value is .25 D-10 meters 


41-56 

D15.8 

TOL3 or R, a nominal value for the 
local truncation error, default value 
is .25 D-7 meters. 


56-70 

D15.8 

STEPMIN, the mmimum stepsize in 
seconds, default value is 5 seconds. 

Note: 

These tolerances must satisfy T| 

< R < T . 
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APPENDIX A2 

ORBGEN PROGRAM DEFINITIONS AND SAMPLE RUN DECK 


Internal Program Definition of Drag Coefficient and Stepsize Option 


ORGBEN provides internally programmed default values and options for-drag 
.coefficient and stepsize when certain orbital conditions occur. 

Drag Coefficient 


C D » the drag coefficient is initialized and position partials will be computed for 
a perigee height less than 800 km. In such a case* 


Cd = 23 


stepsize = 75 seconds 

The input DRAG card may be used to override this option. 
Stepsize 


The nominal stepsize m the fixed step mode is 100 seconds. However, for ec- 
centricity less than .01 and perigee height greater than 1500 km., the stepsize 
is computed as: 

'stepsize = perigee height * 8.849557E-3 + 86.72566 
stepsize maximum = 400 seconds 


The variable Step option is selected for orbits with eccentricity greater than .02. 
In such a case, the initial stepsize is 24 seconds. 

The STEP1, STEP2 and MODE cards of ORBGEN may override these options. 
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SAMPLE RUN DECK SETUP 


//STEP4 EXEC LINKGO, REGION, GO = 400K 

//LINK.SYSLKMDDDSNAME^**** SYSTEM TAPE **** 

//UNK.OBJECT DD * 

ENTRY MAIN 
/* 

//GO.FT35F001DD UNIT = DISK, DCR = (RECFM = VBS.LRECL = 44, BLK 
SIZE = 4404), 

// SPACE = (CYL, (10,1)) 

//GO.DUM DD DSNAME = *.LINK.SYSLMOD, DISP = (OLD,DE LETE) 
//GO.DATA5 DD * 

ORBIT GENERATOR RUN WITH **** CHANGES IN THE STANDARD FORCE 

MODEL AND OUTPUTTING THE PRE- 
DICTED EPHEMERIS EVERY 60 SECS. 


651120 

651126012736.0000 1 
532242.523916446 
-4803.1367878774 
SAT 1.23 

ORBIT 40 651128. 
SOLRAD 1.5 
OUTPUT1 60. 
STEP1 100. 
COEFGP6 


65891 

6357536.01341798 

-3236.3075474982 


172.5 

235500. 


CHG 

C0303 

1.564 

E-07 

CHG 

C1513 

0.0 

E+00 

CHG 

S1513 

0.0 

E+00 

CHG 

Cl 713 

0.0 

E+00 

CHG 

S1713 

0,0 

E+00 


651129 

651129 


END 

EARTH 18183 .986 009000D+14 6378142. 


5275907.36107666 

3657.0878411886 

000500. 

000500. 


298.255 


/* 


// 
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